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The low-energy spectrum of Al-boson clusters with pairwise zero-range interactions is believed to 
be governed by a three-body parameter. We study the ground state of Wboson clusters with infinite 
two-body s-wave scattering length by performing ab initio Monte Carlo simulations. To prevent 
Thomas collapse, different finite-range three-body regulators are used. The energy and structural 
properties for the three-body Hamiltonian with two-body zero-range interactions and three-body 
regulator are in much better agreement with the “ideal zero-range Efimov theory” results than those 
for Hamiltonian with two-body finite-range interactions. For larger clusters we find that the ground 
state energy and structural properties of the Hamiltonian with two-body zero-range interactions 
and finite-range three-body regulators are not universally determined by the three-body parameter, 
i.e., dependences on the specific form of the three-body regulator are observed. For comparison, 
we consider Hamiltonian with two-body van der Waals interactions and no three-body regulator. 

For the interactions considered, the ground state energy of the A-body clusters is—if scaled by the 
three-body ground state energy—fairly universal, i.e., the dependence on the short-range details of 
the two-body van der Waals potentials is small. Our results are compared with the literature. 

PACS numbers: 03.75.-b 


I. INTRODUCTION 

The unitary regime, where the two-body s-wave scat¬ 
tering length is infinitely large, can be reached in ultra 
cold dilute atomic gases using Feshbach resonance tech¬ 
niques [1] . Two-component Fermi gases were realized ex¬ 
perimentally and found to be stable and universal even in 
the large s-wave scattering length regime [2-4], i.e., the 
properties of the system were found to be governed, to a 
very good approximation, by the s-wave scattering length 
as alone and independent of the details of the interac¬ 
tion potential [5-7]. Unitary Bose gases, in contrast, are 
short-lived [8-10]. Their properties depend on the details 
of the interaction potential. Typically, this dependence 
is encapsulated by a three-body parameter [11]. 

Efimov predicted that three identical bosons inter¬ 
acting through two-body potentials with infinitely large 
s-wave scattering length Og and vanishing effective 
range support an infinite number of three-body bound 
states [12]. The binding momenta of the trimers 
{n labels the states) display a geometric scaling, i.e., 
~ 22.6944 [11, 12]. If the binding momen¬ 
tum of one trimer is known, that of the other trimers is 
also known. Importantly, the binding momenta them¬ 
selves cannot be determined solely from a theory that 
is based on two-body zero-range potentials. Rather, a 
three-body parameter is needed to regularize the prob¬ 
lem (i.e., to set the absolute scale of the three-body spec¬ 
trum). The three-body regulator can be introduced in 
many ways. In this work, we consider three different 
regularization approaches: (i) a Hamiltonian with two- 
body zero-range potentials and a zero-range three-body 
potential, (ii) a Hamiltonian with two-body zero-range 
potentials and a purely repulsive three-body potential. 


and (iii) a Hamiltonian with finite-range two-body po¬ 
tentials and no three-body potential. 

Much less is known about four- and higher-body sys¬ 
tems at unitarity [13-21]. 7V-body cluster states are be¬ 
lieved to be attached to each trimer, i.e., for a trimer with 
binding momentum , two A^-body states are believed 
to exist with binding momenta and k^\ 

where and are dimensionless parameters that 
do not depend on n. Whether four- and higher-body 
parameters exist has been under debate in the literature. 

The study of fV-body states attached to Efimov trimers 
is challenging for several reasons. To date, no analytical 
solutions for TV > 4 exist. Numerical treatments have 
to be capable of describing vastly different length scales. 
Eor finite-range two-body interactions, the lowest trimer 
state is typically not a “pure” Efimov state. Thus, one 
would ideally like to investigate fV-body droplets that are 
tied to the first- or second-excited trimer states. The cor¬ 
responding Wbody states {N > 4; see Eig. 1 for an illus¬ 
tration of the four-body spectrum as a function of 1/as) 
are not bound states but resonance states, which are not 
stable with respect to break-up into smaller sub-units. 
Thus, the numerical approach of choice would ideally be 
capable of treating 7V-body resonance states whose size 
is many orders of magnitude larger than the range of the 
underlying two-body potential. 

To bypass these numerical challenges, this work pur¬ 
sues, as have other works before [21, 23], an approach 
that considers fV-body droplets (the thick dashed lines in 
Fig. 1 show the two four-body states) tied to the energet¬ 
ically lowest-lying trimer state (thick solid line in Fig. 1). 
To ensure that the trimer ground state has the key char¬ 
acteristics of a true Efimov trimer state, we employ two- 
body zero-range interactions together with a purely re- 
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FIG. 1: (Color online) Schematic illustration of the en¬ 
ergy spectrum for four identical bosons. The x marks the 
(l/as-^E) = (0,0) point. The dotted line shows the en¬ 
ergy of the weakly-bound dimer. The solid lines show dif¬ 
ferent Efimov trimer states, which become unbound on the 
positive scattering length side at the atom-dimer threshold. 
The dashed lines show “ground state” and “excited state” 
tetramers that are attached to each Efimov trimer. These 
tetramer states hit the dimer-dimer threshold on the positive 
scattering length side (the energy of the two dimers is shown 
by the dash-dotted line). It should be noted that the excited 
tetramer state turns into a virtual state for a certain region 
of positive scattering lengths [22]; this detail is not reflected 
in the plot. 


TABLE I: Summary of potential models considered in this 
work. For each model, the two-body potential V 2 b and the 
three-body potential Vsb are listed. Vib for 2bZR-|-3bZR, 
2bZR-|-3bHC, and 2bZR-|-3bRp is the Fermi-Huang pseu¬ 
dopotential [24]; fls is set to infinity. Vzr(R) for 2bZR+3bZR 
is treated as a zero-range boundary condition. VHC,flo(.R) 
the hardcore repulsive potential; Vhc,Ho(-^) = 0 for R > Rq 
and VHC,i?o(^) = oo for R < Rq. Vq and ro for 2bG, C12 and 
C6 for 2bLJ, CIO and c% for 2bl0-6, and cs and ce for 2b8-6 
are chosen such that the s-wave scattering length is infinitely 
large and the two-body system supports one zero-energy s- 
wave bound state. 


model 


V2b 


Fib 

2bZR-t3bZR 


5(3) (r) 

■^r 

or 

Ezr(R) 

2bZR+3bHC 

^as 

m ° 

5 ( 3 ) (r) 

■S-r 

or 

Lhc,Ro(R) 

2bZR-|-3bRp 

— as 

5 ( 3 ) (r) 

■^r 

dr ^ 

Cp 

RP 

2bG 

Vo exp[ 

-^7(2rg)] 

— 

2bLJ 

C12 _ C6 

j.L'2 


— 

2bl0-6 

CIO _ C6 

lU rpd 


— 

2b8-6 


C6 

■ “ 7^ 


— 


pulsive three-body potential that serves as a regulator; 
we refer to this model as 2bZR-|-3bRp (2b, ZR, 3b, and 
R stand for two-body, zero-range, three-body, and repul¬ 
sive, respectively, and p denotes the power of the repul¬ 
sive three-body potential; see below). The forms of V 2 b 
and Vsb for the model 2bZR-|-3bRp are given in Table I 
and the Hamiltonian H for N particles with mass m and 


position vector rj reads 

N -2 ^ ^ 

^ = -E^^' + E^2b(r,fc)+ E (1) 

j=l j<k j<k<l 

where the two-body potential V 2 b depends on the inter¬ 
particle distance vector rjk {rjk = — r^) and the three- 

body potential Vsb depends on the three-body hyperra¬ 
dius Rjkh 

Rjki = ^(rffc+4+r2;)/3. (2) 

Importantly, the A^-body Hamiltonian H is well behaved, 
i.e., the ground state is well defined thanks to the three- 
body regulator. As we show in Sec. H, the three-body 
regulator produces three-body states that share many 
characteristics with the pure three-body Efimov state. 
Pure three-body Efimov states are obtained if the two- 
body interactions are of zero range and the hyperra- 
dial boundary condition at R 123 = 0 is specified [II]. 
Since the hyperradial boundary condition or logarithmic 
derivative can be imposed via a delta-function in the hy¬ 
perradius, we refer to this model as 2bZR-|-3bZR. 

Our work considers the Af-body ground state using a 
novel Monte Carlo approach [25] that allows for the treat¬ 
ment of two-body zero-range interactions. The Monte 
Carlo approach can unfortunately not treat three-body 
zero-range interactions, i.e., it is not capable of treat¬ 
ing the Hamiltonian 2bZR-|-3bZR. A key objective of the 
present work is then to investigate how the properties of 
A^-body droplets in the ground state, supported by the 
model Hamiltonian 2bZR-|-3bRp, change with the num¬ 
ber of particles and with the power p of the three-body 
regulator. An important question is to which degree the 
A^-body properties are determined by the three-body pa¬ 
rameter. 

For comparison, we also consider Hamiltonian with 
finite-range two-body Gaussian or van der Waals inter¬ 
actions and no three-body interaction. The ground state 
manifolds of these models, referred to as 2bG, 2bLJ, 
2bl0-6, and 2b8-6 (see Table I), lack—as we show—a 
number of key Efimov characteristics. Two-body Gaus¬ 
sian interactions have been employed extensively in the 
literature [19, 21, 26-28], sometimes also in combination 
with a repulsive three-body regulator [23, 29]. 

Although the structural properties of the ground state 
trimers for the Hamiltonian with two-body van der Waals 
interactions differ notably from those for the pure Efi¬ 
mov trimer [30, 31], these systems exhibit universal fea¬ 
tures [27, 32-39]. Specifically, the trimer ground state 
binding momentum at unitarity is, to a good ap¬ 
proximation, determined by the van der Waals length 
Avdw [27, 39] and independent of the short-range de¬ 
tails. For the two-body Lenard-Jones potential, one finds 
Ri 0.230/Lvdw [40], where Tvdw = (v^mc^//i)^/^/ 2 . 
This relationship is nowadays being attributed to van 
der Waals universality. Moreover, the binding momen¬ 
tum spacing of 23.4 between the ground state and the 
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first excited state is quite close to the spacing of 22.6944 
exhibited by consective pure Efimov trimers [40]. It is 
thus interesting to investigate if van der Waals universal¬ 
ity exists for TV > 3, i.e., to answer the question whether 
or not the TV-body ground state energy depends on the 
short-range details of the two-body van der Waals poten¬ 
tial. 

The remainder of this paper is organized as fol¬ 
lows. Section II compares the properties of the three- 
boson system with infinitely large s-wave scattering 
length interacting through 2bZR-|-3bZR, 2bZR-|-3bHC, 
and 2bZR-|-3bRp and illustrates the benefits and limita¬ 
tions of these models. Section III reviews several liter¬ 
ature results for TV-body droplets. Section IV extends 
the calculations for the 2bZR-|-3bRp interaction model 
to clusters with N < 15. In addition to the energy, var¬ 
ious structural properties are discussed in detail. Sec¬ 
tion V compares the results for the model 2bZR-|-3bRp 
with those for systems with two-body hnite-range inter¬ 
actions (i.e., for the models 2bG, 2bLJ, 2bl0-6, and 2b8- 
6 ). Finally, Sec. VI concludes. 


II. THREE-BODY SYSTEM AT UNITARITY 


To understand the three-body system, it is instruc¬ 
tive to rewrite the Hamiltonian H, Eq. (1), for N = 3 
in hyperspherical coordinates [41]. To this end, we first 
separate off the center of mass degrees of freedom and 
restrict ourselves to states with vanishing relative orbital 
angular momentum. For the 2bZR-|-3bZR, 2bZR-|-3bHC, 
and 2bZR-|-3bRp models with infinitely large two-body 
s-wave scattering length a^, the hyperradial and hyper- 
angular degrees of freedom separate [11, 42]. The lowest 
eigen value of the hyperangular Schrodinger equation is 
typically denoted by sq, where sq ~ 1.006i [11, 12]. This 
eigen value enters into the hyperradial Schrodinger equa¬ 
tion with hyperradial Hamiltonian Hu, 


Hr = 




n\sl - 1/4) 

2mi?2 


+ G3b(i?) 


(3) 


(for notational simplicity, the three-body hyperradius 
is denoted by R throughout this section). If V^hiR) 
is equal to zero, the eigen energies of the Hamiltonian 
Hr are not well defined. To make the problem well- 
defined without explicitly introducing a length scale, a 
boundary condition at i? = 0, which serves as a reg¬ 
ulator and defines a scale, can be specified. This is 
the model 2bZR-|-3bZR. The energy spectrum of the 
2bZR-|-3bZR model Hamiltonian displays a perfect ge¬ 
ometric series [11]. For an eigen state with binding mo¬ 
mentum [the corresponding energy is /m], 

there exists a tighter and a looser bound state with bind¬ 
ing momentum Kg” = exp(7r/jso|)K3”^ and Kg"^^^ = 
exp(—Tr/jsoDKg”^ respectively. Here, exp(7r/jso|) is ap¬ 
proximately equal to 22.6944. The three-body spectrum 
for the 2bZR-|-3bZR model is not bounded from below; 
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FIG. 2: (Color online) Breaking of the scale invariance for 
the three-boson system at unitarity with three-body hard¬ 
core regulator. The circles show the difference between the 
binding momentum ratio k'Y '’/of the nth and (n-|- l)th 
states for the model 2bZR-|-3bHC and the ratio exp(7r/|so|) = 
22.6944 for the model 2bZR+3bZR as a function of n. The 
solid line shows a fit to the data points. The breaking of the 
scale invariance becomes weaker with increasing n. 


in our notation, this means that n can take non-positive 
values, i.e., n = ..., —2, —1, 0,1, 2,... There exists an 
infinity of three-body bound states and each hyperradial 
wavefunction YniR) has infinitely many nodes. The hy¬ 
perradial wavefunctions of these states collapse if scaled 
by the binding momentum Kg"^ i.e., R) 

is the same for all states. 

We now consider finite-range three-body regulators. 
As a first toy model, we consider a hardcore repul¬ 
sive three-body potential, i.e., we consider the model 
2bZR-|-3bHC (see Table I). In this case, the hyperan¬ 
gular and hyperradial parts separate as before and the 
Hamiltonian Hr supports a well defined ground state 
with energy or binding momentum Kg^^ (in our no¬ 
tation, n = 1,2,...). For the nth state with binding 
momentum Kg”^, the hyperradial wavefunction has n — 1 
nodes. The circles in Fig. 2 show the difference between 
the binding momentum ratios for the model 2bZR-|-3bHC 
and the model 2bZR-|-3bZR. The binding momentum ra¬ 
tio for the ground and first excited states of the model 
2bZR-|-3bHC is approximately 22.7064. The deviation 
from the model 2bZR-|-3bZR is 0.0120 or 0.053%. As 
we go to excited states, the deviations decrease ex¬ 
ponentially. A log-linear fit of the deviations yields 
Kg^V^s”"*"^^ ~ exp(7 r/jso|) Ri exp(1.823 — 6.244n) (see the 
solid line in Fig. 2). The overlap between the wavefunc¬ 
tion of the ground state of the model 2bZR-|-3bHC and 
the wavefunction of the model 2bZR-|-3bZR with the 
same binding momentum is 0.99947, i.e., the inner re¬ 
gion where the wavefunction for the model 2bZR-|-3bHC 
deviates from the true Efimov wavefunction is insignif¬ 
icant. The three-body hardcore potential breaks the 
scale-invariance and introduces n-dependent energy spac- 
ings. 

The discontinuity of the derivative of the wavefunc¬ 
tion at i? = i?o makes the three-body hardcore regu- 
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FIG. 4: (Color online) Angular distributions for three iden¬ 
tical bosons at unitarity. The circles, triangles and squares 
show the angular distributions Ptot{d), Piniu{d), and Pmax{6) 
for the model 2bZR; these distributions are identical to those 
for the models 2bZR-|-3bHC and 2bZR-|-3bRp. The solid, dot¬ 
ted, and dashed lines show the angular distributions Ptot(0), 
Pmin(^), and Pmax(^) for tile model 2bG. 


FIG. 3: (Color online) Binding momentum characteristics for 
the three-boson system with three-body power law regulator 
at unitarity. The circles show the ratio of the binding mo¬ 
mentum of two consecutive states for the model 2bZR-|-3bRp 
as a function of p. Panel (a) shows the binding momentum 
ratio for the ground and the first excited states while panel 
(b) shows the ratio for the first and the second excited states. 
The solid and dashed lines show the binding momentum ratio 
for the models 2bZR+3bZR and 2bZR+3bHC, respectively. 


lator challenging to treat numerically, at least by the 
path integral Monte Carlo (PIMC) technique employed in 
Sec. IV. Thus, we consider three-body power law poten¬ 
tials, which approach the hardcore potential for p —> oo. 
The circles in Fig. 3 show the binding momentum ratios 
for the model 2bZR-|-3bRp as a function of p. Figures 
3(a) and 3(b) show the binding momentum ratios for the 
ground and first excited states, and the first and sec¬ 
ond excited states, respectively. As expected, the bind¬ 
ing momentum ratios approach the value for the model 
2bZR-|-3bHC (dashed lines) in the large p limit. For com¬ 
parison, the solid lines show the binding momentum ratio 
for the model 2bZR-|-3bZR. The deviations between the 
binding momentum ratios for the 2bZR-|-3bRp and the 
2bZR-|-3bHC models are largest for p = 4. Similar to the 
model 2bZR-|-3bHC, the binding momentum ratios for 
the model 2bZR-|-3bRp approach the value exp(7r/|so|) 
exponentially with increasing n. 

The spacing of the momenta is not the only way to 
characterize how universal the system is, i.e., how close a 
given system is to the true Ehmov scenario described by 
the model 2bZR-t-3bZR. The structural properties pro¬ 
vide additional insights. Indeed, the structures of weakly- 
bound three-body systems with positive have recently 
been measured [30, 31]. We first look at the distribu¬ 
tion of the angles Ojki between each pair of position vec¬ 
tors, 6jki = arccos(fjfc • ffcj). The distribution Ptot{0) 


considers all three angles of each triangle, while the dis¬ 
tribution Pmin(^*) [Tmax(6*)] considers only the smallest 
[largest] of the three angles of each triangle. The nor¬ 
malizations are chosen such that Ptot (9)d9 = 3 and 
Jq Pmin{9)d9 = /J" Pmax(6')d6» = I. For infinitely large 
Os (as considered throughout this section), these angular 
distributions only depend on the hyperangles and not on 
the hyperradius. Thus, they are the same for the mod¬ 
els 2bZR-b3bZR, 2bZR-b3bHC, and 2bZR-b3bRp. The 
circles, triangles, and squares in Fig. 4 show Ptot(^), 
Pmin{9), and Pina,x{9), respectively, for these models. 
Ptot{9) is approximately linear and approaches a finite 
value for 0 —>■ 0. We are interested in the angular distri¬ 
butions for two reasons, (i) For the models 2bG, 2bLJ, 
2bl0-6, and 2b8-6, the hyperangular and hyperradial de¬ 
grees of freedom do not separate and the difference be¬ 
tween their angular distributions and those for the two- 
body zero-range models provides valuable insights (see 
Ref. [40]). (ii) For the iV-body clusters, the angular dis¬ 
tributions, which depend on both the hyperangles and 
the fV-particle hyperradius, can serve to monitor the 
three-body correlations. 

Solid, dotted, and dashed lines in Fig. 4 show the an¬ 
gular distributions Ptat[9), T’min(^), and Pmax(^*), respec¬ 
tively, of the three-body ground state for the model 2bG. 
Compared to that for the two-body zero-range models, 
the angular distribution near 0 = 0 for the finite-range 
model displays distinctly different behavior. For the 
Gaussian model, the probability of Ending an angle of 
zero is zero and the angular distribution peaks at around 
O.ITtt or 31°. For the zero-range model, the angular dis¬ 
tribution peaks at 0 and Ptot(O) is finite. This is because 
the zero-range boundary condition makes the probabil¬ 
ity to find two particles at the same position finite. A 
vanishing interparticle distance corresponds to a triangle 
in which one of the three angles 9jki is zero. Since the 
angular distributions for the models 2bZR-|-3bZR and 
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FIG. 5: (Color online) Radial density p(r) for three identical 
bosons at unitarity (r is measured relative to the center of 
mass of the three-body system). The dashed and solid lines 
show p(r) for the models 2bZR-|-3bZR and 2bZR-|-3bRp with 
p — 6, respectively. 


2bG show distinctly different features, one might expect 
that the binding momentum ratios for these two 

models also differ. The value of for the model 

2bG is 22.983, which differs by only 1.27% from that for 
the model 2bZR-|-3bZR. This indicates that it is insuffi¬ 
cient to only evaluate the binding momentum ratios to 
judge how universal the system is. We note that the dis¬ 
tribution P{6) for the ground state of the TV = 3 system 
with two-body Lenard-Jones interactions is quite similar 
to that for the ground state of the TV = 3 system with 
two-body Gaussian interactions [40]. 

We now consider the radial density p(r) [r is mea¬ 
sured relative to the center of mass of the three-body 
system) for the models 2bZR-|-3bZR and 2bZR-|-3bRp 
with p = 6. The radial density p(r) is normalized such 
that 47r p{r)r‘^dr = TV and depends on the hyperra¬ 
dius and the hyperangles. The dashed and solid lines 
in Fig. 5 show the radial density p{r) for the models 
2bZR-|-3bZR and 2bZR-|-3bRp with p = 6, respectively. 
For the latter, the ground state density is shown. The 
radial densities are scaled by their respective binding mo¬ 
mentum ^ 3 . The solid and dashed lines agree well in the 
large r region and differ notably in the small r region. 
The deviation in the small r region comes from the fact 
that the hyperradial density for the model 2bZR+3bZR 
decays much slower for small R than that for the model 
2bZR-|-3bRp. Note that even though the radial densities 
for the two models differ by about a factor of two in the 
small r region, the difference between the integrated con¬ 
tributions is small because the volume element contains 
an factor. 



FIG. 6: (Golor online) Energy per particle of TV-boson clusters 
at unitarity. (a) Summary of literature results. The dashed 
and dotted lines show the analytical prediction by Gattobi- 
gio and Kievsky [21] and Nicholson [43], respectively. The 
triangles show the diffusion Monte Carlo (DMC) energies for 
a Hamiltonian with two-body square well interaction and re¬ 
pulsive three-body hardcore regulator [23]. The diamonds 
show the energy for the model 2bG [28]. (b) Summary of our 
PIMC calculations. The circles and pluses are for the model 
2bZR-|-3bRp with p = 4 and 8, respectively; the error bars 
(not shown) are of the order of the symbol sizes. The squares, 
diamonds, and triangles are for the model 2bZR+3bRp with 
p = 5, 6, and 7, respectively; the error bars (not shown) are 
smaller than the symbol sizes, (c) Summary of our calcula¬ 
tions for two-body van der Waals models. The circles, crosses, 
and squares show our DMC results for the models 2bLJ, 2bl0- 
6, and 2b8-6, respectively. 


III. N-BODY CLUSTERS AT UNITARITY: 
OVERVIEW OF LITERATURE RESULTS 

This section discusses various literature results for the 
energy of weakly-bound TV-body droplets (TV > 3) at 


unitarity. The diamonds in Fig. 6(a) show the TV-boson 
energy per particle E^/N for the model 2bG as a func¬ 
tion of TV [21, 26, 28]. The energy per particle increases 
approximately linearly with TV for TV > 6 (for smaller TV, 
some deviations from the linear behavior exist). Based 
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on the fact that the energy per particle, and correspond¬ 
ingly the binding momentum, scale linearly with N for 
the model 2bG, Gattobigio et al. [21] proposed an ana¬ 
lytical form for the A^-boson system with two-body zero- 
range interactions and fixed three-body parameter, 

^ = (4) 

[see the dashed line in Fig. 6 (a)]. The ratio K 4 /K 3 is not 
taken from the ground state calculations for the Gaussian 
two-body interaction model, for which K 4 /K 3 = V5.86, 
but from Deltuva’s calculations for highly excited four- 
body resonance states. Deltuva finds the universal ratio 
kaIK 3 = V4.61 [17]. Gattobigio et a/.’s expression, con¬ 
verted to the energy, exhibits a leading order and 
sub-leading order N dependence. 

It should be noted that the ground state energy of the 
Hamiltonian with pairwise Gaussian interactions scales 
differently with for iV > 10 than that of Hamilto¬ 
nian with pairwise interactions with short-range repul¬ 
sion. For interactions with repulsive core, it is well- 
established that the energy per particle increases weaker 
than linear for A^ > 10 (see, e.g., the literature on he¬ 
lium and tritium droplets [44, 45]). Gattobigio et al. [26] 
noted that Eq. (4) applies not only to systems with zero- 
range interactions but also to systems with finite-range 
interactions in the regime where E/N is approximately 
proportional to N (e.g., to helium droplets with N < 10). 
In this case, the ratio for the hnite-range potential 

is taken as input and the binding momentum for A'^ > 4 
is predicted. We return to this discussion in Sec. V. 

Independent evidence for the leading-order N depen¬ 
dence of the energy per particle for the Hamiltonian with 
two-body zero-range interactions comes from lattice cal¬ 
culations for even N [43] . Assuming that the distribution 
of the two-body correlator is exactly log normal, Nichol¬ 
son deduced an analytical expression for the energy per 
particle, Em/N = {N/2 — l)Ei/^ [see the dotted line in 
Fig. 6 (a)] [43]. To plot this expression, we used Deltuva’s 
value of A 4 /A 3 = 4.61. It should be noted that the coef¬ 
ficients predicted by Gattobigio et al. and Nicholson for 
the leading order N dependence differ by about a factor 
of 2 . 

A somewhat different A-dependence of the energy per 
particle was observed in the numerical calculations by 
von Stecher [see triangles in Fig. 6 (a)] [23]. In fact, 
the idea to use a three-body regulator, as in our model 
2bZR-|-3bRp, to make the ground state trimer large and 
Efimov-like was introduced in Ref. [23]. Von Stecher em¬ 
ployed a model Hamiltonian with two-body square well 
potential with infinitely large two-body s-wave scattering 
length and three-body hardcore potential. For N < 10, 
the energy per particle increases approximately linearly 
with increasing N. For larger A, the triangles in Fig. 6 (a) 
flatten. Reference [46] interpreted this as a turnover to 
a A° dependence of the energy per particle. Such a be¬ 
havior suggests a saturation of the density for large A. 
This saturation would be a consequence of the balance 


of the two-body attractive and three-body repulsive in¬ 
teractions. 

The discussion above shows that the dependence of 
the energies tied to Efimov trimers is not well under¬ 
stood. Specifically, neither the functional form of the en¬ 
ergy per particle nor the coefficients are agreed upon. In 
the following sections, we attempt to understand where 
the discrepancies of the literature results come from. 


IV. A-BODY RESULTS AT UNITARITY FOR 
THE MODEL 2bZR+3bRp 

To calculate the A-boson energy for the Hamiltonian 
with interaction model 2bZR-|-3bRp, we apply the PIMC 
technique [25, 47]. The PIMG technique is an, in prin¬ 
ciple, exact finite-temperature method; the errors, which 
originate from the discretization of the imaginary time 
and the stochastic evaluation of integrals, can be re¬ 
duced systematically. To obtain the ground state en¬ 
ergy of the A-boson Hamiltonian, the PIMC approach 
has to be extended to the zero-temperature limit. Typ¬ 
ically, this is achieved by the path integral ground state 
approach [47, 48]. Here, we pursue an alternative strat¬ 
egy. Namely, we work in the finite temperature regime 
where the thermal contribution to the energy is known 
and where the structural properties of interest are not 
affected by the temperature. This approach was intro¬ 
duced and benchmarked in Ref. [28]. The basic idea is 
to place the droplet in a weak external harmonic confine¬ 
ment, whose angular frequency to is chosen such that the 
center of mass energy spectrum becomes discretized and 
the relative motion is unaffected by the trap. This re¬ 
quirement corresponds to \Em\ ^ fno. Since the density 
of states of the harmonically trapped center of mass pseu¬ 
doparticle is known analytically, the ground state energy 
Em of the A-boson droplet in free space can be extracted 
from the finite-temperature energy [25, 28]. 

The circles, squares, diamonds, triangles, and pluses 
in Fig. 6 (b) show the energy per particle for the model 
2bZR-|-3bRp with p = 4, 5, 6 , 7, and 8 , respectively, as 
a function of A (see also Table H and the Supplemen¬ 
tal Material [54]). For each p, the energy per particle is 
scaled by the respective trimer energy per particle. For 
a fixed p, the energy per particle increases monotonically 
and smoothly as a function of A, i.e, even-odd effects, 
which have been observed in trapped and homogeneous 
two-component Fermi gases [49-51], are—if existent— 
smaller than our statistical error bars. For fixed A, 
the scaled energy per particle increases with increasing p 
(p > 4); this increase becomes smaller with increasing p. 
Similarly to von Stecher’s energy per particle [23] [trian¬ 
gles in Fig. 6 (a)], the scaled energy per particle increases 
roughly linearly for smallish A and then flattens out for 
larger A. This effect is most pronounced for p = 4 and 
5, where the flattening sets in around A = 8 — 10, and 
least pronounced for p = 8 . The reason for the flattening 
is that the clusters develop, for sufficiently large A, more 
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TABLE II: PIMC energies for the model 2bZR+3bRp 
for N = 4 — 15. Columns 2-4 show the scaled energy 
En/N/{E 3 / 3 ) for p = 5, 6, and 7, respectively. The error 
bars (not explicitly reported) are around 3%. 

N 2bZR+3bR5 2bZR-b3bR6 2bZR-b3bR7 


4 

3.46 

3.64 

3.73 

5 

6.19 

6.53 

6.70 

6 

8.69 

9.42 

9.81 

7 

10.9 

12.0 

12.6 

8 

12.8 

14.3 

15.1 

9 

14.5 

16.4 

17.5 

10 

15.9 

18.3 

19.7 

11 

17.3 

20.0 

21.5 

12 

18.4 

21.5 

23.3 

13 

19.4 

22.8 

25.0 

14 

20.3 

24.2 

26.4 

15 

21.1 

25.2 

27.8 



P 

FIG. 7: (Color online) Comparison of our PIMC energies 
(left) and literature results (right) for selected N. Panels (a), 
(b), and (c) show our PIMC energy per particle for 77-boson 
clusters interacting through the model 2bZR+3bRp as a func¬ 
tion of p for N — 6, 10, and 13, respectively. For comparison, 
panels (d), (e), and (f) show the energy per particle from 
the literature for the same N. The triangles, diamonds, and 
squares are from von Stecher [23], Nicholson [43], and Gatto- 
bigio et al. [21]. Since the work by Nicholson is restricted to 
even N, comparison for = 13 cannot be made. 

than one pair distance scale (see below for more details). 

The circles in Fig. 7 replot the PIMC energy per par¬ 
ticle for selected N. As the power p increases, the scaled 
energy approaches a constant. Based on our discussion 
in Sec. II, the p —>■ 00 energy should coincide with the en¬ 
ergy for the model 2bZR-|-3bHC. It is thus instructive to 
compare our scaled energies, extrapolated by eye to the 
p —>■ 00 limit, with those obtained by von Stecher [23], 
who employed a two-body square well potential and a 
three-body hardcore regulator [see triangles in Figs. 7(d)- 


7(f)]. We find that our p —>■ 00 energy per particle lies 
above von Stecher’s energy per particle by something like 
10 - 20%, 20 - 30%, and 30 - 50% for N = 6, 10, and 13, 
respectively. Since the three-body sectors are treated on 
consistent footing (3bRp—j-dbHC as p — 00 ), we specu¬ 
late that the difference arises from the different two-body 
interactions. However, we did not perform calculations 
to conhrm this and can thus not rule out other reasons. 
As can be seen from Fig. 7, Nicholson’s energy prediction 
lies notably below our large p energies while Gattobigio et 
al.’s prediction lies above our large p energies for N >8. 

If the iV-body energies were determined solely by a 
three-body parameter K 3 , the model 2bZR-|-3bRp for dif¬ 
ferent p would yield the same scaled energies, i.e., the 
symbols in Fig. 6 (b) would collapse to a single curve. The 
fact that they do not collapse indicates that the three- 
body parameter is not sufficient to predict the energy of 
the A^-boson clusters, at least not for the models consid¬ 
ered. To gain more insight into this, it is instructive to 
analyze the length scales of the model 2bZR-|-3bRp. Four 
length scales can be identified (see rows 3-6 of Table III), 
(i) The characteristic length scale Lp of the three-body 
repulsive potential, (ii) The length scale Z 3 defined by 
the three-body binding energy, (iii) The length scale Ln 
defined by the energy of the cluster. And, (iv) the length 
scale In associated with the energy per particle of the 
cluster. Inspection of the definitions given in Table III 
shows that Ln and In are not independent. 

Forp = 4-8, we find L^/Lp r; 29.3,28.8,27.6,26.6, and 
25.9, i.e., the trimer is significantly larger than the scale 
of the underlying repulsive three-body potential. This 
ensures, as discussed in Sec. II, that the trimer ground 
state described by the model 2bZR-|-3bRp with p > 4 
exhibits the key characteristics of an Efimov state. It is 
instructive to alternatively think about the trimer size in 
terms of the average interparticle spacing f. For trimers 
with p = 4-8, we find f/Lp « 18.7,18.5,17.7, 17.1, and 
16.6. 

For p = 6 , we find that Ln/Lp changes from 11.2 for 

= 4 to 8.37 for AT = 5 to 2.46 for N = 15. This sug¬ 
gests that the A^-boson droplet “sees” increasingly more 
of the three-body regulator as N increases, i.e., that the 
dependence of En/N on p increases with increasing N. 
The length scale In, in contrast, suggests a larger sepa¬ 
ration of scales; for N = 13, e.g., we have In/Lp = 7.69 
for p = 4 and In/Lp = 6.85 for p = 8 . In fact, if En/N 
scales as N, then Ln and In scale as 1/A^ and I/'/N, 
respectively. If En/N scales as N^, then Ln and In 
scale as l/^/lL and N^, respectively. This implies that— 
unless the energy scales linearly (or even weaker) with 
N for large N —the properties of the Al-boson droplets 
are expected to be notably affected by the choice of the 
three-body regulator. 

Alternatively, one can consider the average interparti¬ 
cle distance r and the average sub-three-body hyperra¬ 
dius R. The squares in Fig. 8 (a) show the average inter¬ 
particle spacing r, i.e., the expectation value of the pair 
distance, as a function of N in units of I/K 3 (left axis) 














TABLE III: Summary of the definitions of length scales. The van der Waals length Lvdw is defined in Ref. [1]. Lp for p — 6 
agrees with Lvdw if rn is replaced by the reduced two-body mass m/ 2 . 


length scale 

definition 

description 

L, 

ro 

characteristic length scale of the two-body Gaussian potential 

LvdW 

{^mcz/hY ^'^/2 

characteristic length scale of the two-body van der Waals potential 

Lp 

[l/(p-2)y2^/fi]2/(P-2) 

characteristic length scale of the three-body repulsive potential 

L 3 

1/k3 = h/^Jm\Ez\ 

length scale set by the three-body binding energy 

Ln 

l/«:jv = h/^m\EN\ 

length scale set by the A-body binding energy 

In 

h/^Jm\EN\/N = y/NLN 

length scale set by the A-body binding energy per particle 

f 


average interparticle spacing 

R 


average sub-three-body hyperradius 



FIG. 8 : (Color online) Expectation value f of the pair dis¬ 
tance as a function of N for A^-boson systems interacting 
through various models, (a) The squares are for the model 
2bZR+3bRp with p = 6 . (b) The triangles are for the model 
2bG. (c) The circles are for the model 2bLJ. The error bars 
show the variance of the pair distance. The pair distances 
are plotted using two different units: (i) the inverse three- 
body binding momentum (left axis) and (ii) the characteristic 
length scale of the model Hamiltonian (right axis). 



FIG. 9: (Color online) (a) Expectation value R of the sub- 
three-body hyperradius (triple size) as a function of N for N- 
boson systems interacting through the model 2bZR-|-3bR6. 
The error bars show the variance of the triple size, (b) 
Triple distribution function Ptripie(R) for the = 13 clus¬ 
ter scaled using the three-body binding momentum kz- The 
solid lines from top to bottom at k,zR = 0.6 are for the model 
2bZR+3bRp with p = 4, 5, 6 , 7, and 8 . The inset replots the 
triple distribution functions using the binding momentum K 13 
of the A = 13 droplet. In these units, the triple distribution 
functions for different p collapse. 


and in units of (right axis) for the model 2bZR-|-3bRp 
with p = 6. The error bars indicate the variance Ar of 
the pair distance, Ar = \/ {r"^) — where () indicates 
the quantum mechanical expectation value [55]. As the 
number of particles N increases, both the mean and vari¬ 
ance of the pair distance are nearly constant. The mean 


and variance of the pair distance are about one order of 
magnitude larger than the internal length scale Lp. The 
relatively large variance of the Hamiltonian with model 
interaction 2bZR-|-3bRp implies that the clusters are dif¬ 
fuse and liquid-like. The squares in Fig. 9(a) show the av¬ 
erage sub-three-body hyperradius R, i.e., the expectation 
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FIG. 10: (Color online) (a) Maximum density pmax as a func¬ 
tion of N for A^-boson systems interacting through various 
models. The circles, squares, and diamonds are for the model 
2bZR-|-3bRp with p = 5 (lowest data set), 6, and 7 (highest 
data set), respectively. For comparison, the line is for the 
model 2bG. (b) Same data as in (a) but replotted as the min¬ 
imum average interparticle distance (pmax)”^^®. The right 
axis shows the data for the model 2bZR-|-3bR6 in units of Le. 


value of the triple size, as a function of N for the model 
2bZR-|-3bRp with p = 6 . The error bars indicate the 
variance. The mean and variance of the sub-three-body 
hyperradius behave similar to the mean and variance of 
the pair distance. 

The average pair distance and sub-three-body hyper¬ 
radius are obtained by averaging over all possible pairs 
and triples regardless of whether or not the particles are 
close to each other. To get more “local” information, 
we calculate the maximum density and subsequently the 
closest pair distance. The circles, squares, and diamonds 
in Fig. 10(a) show the maximum pmax of the radial den¬ 
sity for the model 2bZR-t-3bRp with p = 5,6, and 7, 
respectively, as a function of N. We find that unlike for 
N = 3 (see Fig. 5), the radial density peaks at r = 0 
for iV > 4. For all p, the maximum of the radial den¬ 
sity is roughly a constant for the largest N considered. 
This constant depends—as the energy per particle—on 
the three-body regulator. The circles, squares, and dia¬ 
monds in Fig. 10(b) show the smallest average pair dis¬ 
tance for the model 2bZR-|-3bRp with p = 5, 6 , and 7, 
respectively, as a function of N. The smallest average 
pair distance decreases with increasing N and approxi¬ 
mately saturates for the largest N considered. The small¬ 
est average pair distance is only about five times larger 
than the characteristic length scale Lp of the three-body 
regulator. 

The above length scale discussion can be expanded by 
considering distribution functions. The scaled pair dis¬ 
tribution function 47 rr^Ppair(r’), normalized according to 



FIG. 11: (Color online) Scaled pair distribution function 
47 rr^Ppair(r) for 77 = 13 bosons interacting through various 
models, (a) The solid lines from top to bottom at rear = 0.8 
are for the model 2bZR-|-3bRp with p — 4-8, scaled using 
the three-body binding momentum res. The inset replots the 
pair distribution functions scaled using the binding momen¬ 
tum rei 3 of the 77 = 13 droplet. In these units, the pair 
distribution functions for different p collapse, (b) The dashed 
and dotted lines show the scaled distribution functions for 
the models 2bLJ and 2bG, respectively, using the three-body 
binding momentum res. 


/o°°= 1, tells One the probability to find 
two particles at a distance r from each other. The lines 
from top to bottom at = 0.8 in Fig. 11(a) show the 
scaled pair distribution function 47rr^Ppair('r) for 77 = 13 
interacting through 2bZR-|-3bRp with p = 4-8. The am¬ 
plitude at r = 0 is finite and roughly independent of p. 
This makes sense as it is a signature of the two-body 
zero-range interactions, which enforce a finite amplitude 
at r = 0 . 

The triple distribution function Ptripie(fi), normalized 
according to Ptripie{R)dR = 1 , tells one the probabil¬ 
ity to find three particles with sub-three-body hyperra¬ 
dius R. The solid lines from top to bottom at K 3 R = 
0.6 in Fig. 9(b) show the triple distribution function 
-Ptripie(.R) for 77 = 13 interacting through 2bZR-|-3bRp 
with p = 4-8. The triple distribution functions are broad 
and structureless, indicating that the clusters are diffuse 
and liquid-like and that no small three-body sub-systems 
are formed. 

Figures 9(b) and 11(a) show that the distribution func¬ 
tions Ppaii{r) and Ptripie{R) do not collapse if scaled by 
the three-body binding momentum K 3 . The distribution 
functions for p = 4 are notably broader than those for 
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FIG. 12: (Color online) Angular distribution Ptot{d) for N- 
boson clusters interacting through the model 2bZR+3bRp 
with p = 6 . The lines from top to bottom at 0 = 0 are 
for N = 5,6, 7, 9, and 13. 

p > A. Figures 9(a) and 11(a) suggest that the distribu¬ 
tion functions converge in the large p limit (i.e., in the 
three-body hardcore regulator limit). Similar behavior is 
observed for other N. As shown in the insets of Figs. 9(b) 
and 11 (a), the distribution functions collapse to a very 
good approximation to a single curve if scaled by the 
binding momentum kn of the A^-body droplet. This can 
be understood as a new type of universality, which is 
weaker than the “Efimov universality”: The binding mo¬ 
mentum KN allows one to collapse the distribution func¬ 
tions for the models 2bZR-|-3bRp for sufficiently large p 
but kn is not determined by K 3 (the latter would consti¬ 
tute “Efimov universality”). The dominance of /cat arises 
because the vast majority of the wave function amplitude 
is located in the classically forbidden region [52] (for pure 
zero-range interactions, the classically allowed region is 
reduced to a single point). 

At the three-body level, the angular distribution func¬ 
tions for the models 2bZR-|-3bRp and 2bZR-|-3bZR co¬ 
incide since the hyperradial and hyperangular degrees of 
freedom separate. This is not the case for N > 3, since 
the three-body regulator depends on the 7V-body hyper¬ 
radius and a subset of the 37V — 4 hyperangles. For fixed 
N, we find that the dependence of the angular distribu¬ 
tion functions Ptot(^) on the power p of the three-body 
regulator is small [much smaller than the dependence of 
Ppahir) and /triple(7?) on p]. Figure 12 shows the an¬ 
gular distribution function Ptot(^*) for iV-boson clusters 
interacting through 2bZR-|-3bR6 for various N. The lines 
from top to bottom at 0 = 0 are for N = 5,6, 7,9, and 
13. As the number of particles increases, the probability 
of finding triangles with small angles decreases but re¬ 
mains finite. Intuitively, this is because /tot(^) accounts 
for all the trimer configurations and not just the “closest 
trimers”. 

Combining the information displayed in Figs. 6-12, 
the key characteristics of the ground state of A^-boson 
droplets interacting through the model 2bZR-|-3bRp with 
p > A can be summarized as follows: (i) The depen¬ 
dence of the energy and the structural properties on the 


three-body regulator decreases with increasing p; (ii) the 
dependence of the energy and the structural properties 
on the three-body regulator cannot be explained by sim¬ 
ple length scale arguments (the separation of scales is 
largest for the p = A regulator and smallest for the p = 8 
regulator); (iii) the pair and triple distribution functions 
collapse to a very good approximation to a single curve 
if scaled by the binding momentum of the A^-body sys¬ 
tem, suggesting that 1/kn and not I/K 3 is the governing 
length scale for A^ > 3. 


V. RESULTS FOR OTHER INTERACTION 
MODELS 

We now compare the findings for A-boson systems in¬ 
teracting through the model 2bZR-|-3bRp with p = 4 — 8 
(see the previous section) with those for A-boson sys¬ 
tems interacting through the models 2bG, 2bLJ, 2bl0-6 
and 2 b 8 - 6 . 

We start our discussion with the model 2bG, for which 
the energy per particle scales, to a very good approx¬ 
imation, linearly with N for N > 6 [see diamonds in 
Fig. 6 (a)]. The model 2bG has no repulsive core and is 
characterized by a single length scale, the width tq. Us¬ 
ing a simple variational Gaussian product wave function 
in the single-particle coordinates, one can readily show 
that the ground state energy scales as and that the 
peak density increases quadratically with N. Indeed, our 
calculations shown in Figs. 8 (b) and 10 for up to A^ = 15 
clearly support that the droplet shrinks with increasing 
N. As can be seen in Fig. 8 (b), the average interpar¬ 
ticle distance quickly decreases to a value smaller than 
tq. We conclude that the N'^ scaling of the energy for 
the model 2bG predominantly reflects the absence of a 
repulsive core in the potential energy and less so Efimov 
characteristics. 

Next, we discuss the properties of the Hamiltonian in¬ 
teracting through the van der Waals models 2bLJ, 2bl0- 
6 , and 2b8-6. Our calculations at unitarity are per¬ 
formed using the same atomic mass and the same Cg 
coefficient for the three models while the short-range 
coefficients are tuned such that the dimer supports a 
single s-wave bound state with zero energy. For the 
three-body system, we find KsTvdw = 0.230 for the 
model 2bLJ, K 3 Lvdw = 0.233 for the model 2bl0-6, and 
KsUvdw = 0.245 for the model 2b8-6, i.e., the three-body 
binding momentum depends weakly on the short-range 
scale of the two-body potential. The A-body energies per 
particle, in units of the three-body energy per particle, 
are summarized in Table IV. These energies are obtained 
by the DMG approach [53]. Dividing the V-body ener¬ 
gies by the corresponding three-body energy, the energy 
per particle curves for the three van der Waals interac¬ 
tion models nearly collapse [see Fig. 6 (c)]. This can be 
interpreted as van der Waals universality in the V-body 
sector. Due to the repulsive core, the energy per particle 
flattens around N = 10, indicating that the system starts 
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TABLE IV: DMC energies for the Hamiltonian with two- 
body van der Waals interactions for V = 4 — 15. Columns 2-4 
show the scaled energy En/N /for the models 2bLJ, 
2bl0-6, and 2b8-6, respectively. The error bars (not explicitly 
reported) are around 1%. 


A 

2bLJ 

2bl0-6 

2b8-6 

4 

3.978 

3.953 

3.960 

5 

7.827 

7.841 

7.887 

6 

11.95 

11.99 

12.12 

7 

16.07 

16.15 

16.40 

8 

20.09 

20.24 

20.59 

9 

23.94 

24.15 

24.69 

10 

27.57 

27.89 

28.57 

11 

31.07 

31.44 

32.29 

12 

34.37 

34.81 

35.86 

13 

37.50 

38.02 

39.25 

14 

40.46 

41.06 

42.41 

15 

43.27 

43.97 

45.46 


0.1 


D. 

V - 0.1 
- 0.2 

04 8 12 16 

N 

FIG. 13: (Color online) Assessing the applicability of Eq. (4) 
for A-boson systems with two-body finite-range interactions 
at unitarity. Circles and triangles show the normalized differ¬ 
ence — kn)/kn for the models 2bG and 2bLJ, respec¬ 

tively, as a function of the number of particles N. 



VI. CONCLUSIONS 


to grow outward, i.e., starts to form a “second layer” (of 
course, the system is liquid-like and individual layers can¬ 
not be identified). Consistent with this. Fig. 8 (c) shows 
that the average interparticle distance first decreases with 
increasing N and then slowly increases for V > 8 . 

The dashed line in Fig. 11(b) shows the pair distribu¬ 
tion function of the N = 13 system interacting through 
the model 2bLJ. The amplitude in the small r region 
is suppressed compared to the other interaction mod¬ 
els considered due to the repulsive two-body core. Scal¬ 
ing r^Ppair(c) using K 13 (not shown) does not bring the 
pair distribution function for the model 2bLJ in agree¬ 
ment with the scaled pair distribution functions shown 
in the inset of Fig. 11(a) for the model 2bZR-|-3bRp with 
p = 4 — 8 . This reflects the fact that a notably smaller 
fraction of the wave function amplitude resides in the 
classically forbidden region for the model 2bLJ than for 
the model 2bZR-|-3bRp with p = 4 — 8 . 

As already mentioned in Sec. Ill, Eq. (4) applies, ac¬ 
cording to Ref. [26] , not only to systems with zero-range 
interactions but also to systems with finite-range two- 
body interactions. To assess the applicability of Eq. (4), 
we denote the left hand side of Eq. (4) by jK 3 

and plot the normalized difference between /K 3 

and the exact Kff /K 3 , as determined by our calcula¬ 
tions. Circles and triangles in Eig. 13 shows the quantity 
(k^^” — kn)/kn for the models 2bG and 2bLJ, respec¬ 
tively. Eor N = 3 and A = 4, the normalized difference 
is zero by construction. For A > 4, the normalized dif¬ 
ference is negative for the model 2bG and positive for 
the model 2bLJ. The deviations from the functional form 
proposed by Gattobigio et al. increase roughly linearly 
with A for the model 2bLJ, reaching 13% for A = 15, 
and non-linearly for the model 2bG, reaching —20% for 
A = 15. Thus if high accuracy predictions are sought, 
then Eq. (4) should be used with caution. 


This paper studied weakly-bound Bose droplets at uni¬ 
tarity. These systems are obtained by adding one atom 
at a time to an Efimov trimer or a weakly-bound trimer 
with Efimov characteristics. We carefully analyzed the 
three-body system and then studied larger systems. 

The three-body ground state of the Hamiltonian with 
two-body zero-range interactions and repulsive three- 
body potential (model 2bZR-|-3bRp) is a nearly ideal Efi¬ 
mov state. The premise was (see also Ref. [23]) that this 
would allow us to determine the universal properties of 
droplets tied to a three-body Efimov state by studying 
A-body ground states. Somewhat surprisingly, we found 
dependences of the ground state cluster properties on the 
three-body regulator, suggesting that the ground states 
become less universal with increasing A. This is a some¬ 
what disappointing finding as the treatment of A-body 
excited and resonance states, which are expected to ex¬ 
hibit universal characteristics, is a computationally much 
more demanding task. Yet, our study revealed a dif¬ 
ferent type of universality for these model Hamiltonian. 
We found that if the lengths are scaled by the A-body 
binding momentum, then the dependence on the three- 
body regulator diminishes notably. This suggests that 
the ground states of these systems are halo states [52], 
i.e., states whose amplitude is predominantly located in 
the classically forbidden region. The A-body binding 
momentum itself is, however, not—as it would be in the 
case of A-body Efimov universality—determined by the 
three-body binding momentum, especially not as A in¬ 
creases. 

Hamiltonian with two-body van der Waals interaction 
at unitarity were also investigated. It was found that 
the energy per particle, if scaled by the three-body en¬ 
ergy, collapses to a very good approximation to a sin¬ 
gle curve, suggesting that the short-range details of the 
van der Waals interaction impact the three- and higher- 
body sectors in a similar manner (i.e., the short-range 
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details are to a very good approxiniation “taken out” 
by scaling by the three-body energy). The calculations 
presented were for Lenard-Jones and modified Lenard- 
Jones potentials; the latter potentials have a —ce/r® tail 
but a softer repulsive core at small distances than typical 
van der Waals interactions. We also performed calcula¬ 
tions for (i) the true helium-helium potential scaled by 
an overall factor such that the s-wave scattering length is 
infinitely large and (ii) the true helium-helium potential 
with modified short-range potential such that the s-wave 
scattering length is infinitely large (these models were la¬ 
beled He-He(scale) and He-He(arctan) in Ref. [40]). The 
energy per particle curves for these systems, which have 
a more complicated long-range tail, also collapse, to a 
very good approximation, to the same curves as those 
for 2bLJ, 2bl0-6, and 2b8-6 if scaled by the three-body 
energy. The structural properties, specifically the pair 
and triple distribution functions, for the van der Waals 
systems do not collapse to the same curves as those for 
the 2bZR-|-3bRp model with p = 4 — 8 if scaled using 
the fV-body binding momentum kat, suggesting that a 
good portion of the wave function amplitude of the van 
der Waals systems is located in the classically allowed 


region. 

In the future, it would be interesting to extend the cal¬ 
culations presented here to excited and resonance states. 
We expect that the A^-body properties become univer¬ 
sal if sufficiently high excitations are being considered. 
In the four-body sector, e.g., Deltuva [17] extracted the 
universal numbers for K 4 /K 3 by going to high-lying res¬ 
onance states (in this case, “high-lying” means third or 
higher resonance states). Extending calculations such as 
those conducted by Deltuva to TV > 4 is, however, chal¬ 
lenging. It would also be interesting to extend the studies 
presented in this paper to finite s-wave scattering lengths 
and to Bose droplets with an impurity. 
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The PIMC energies reported in Table II of the main text were obtained using the “fourth-order propagator” 
described in Ref. [I] using a fixed time step of 0.00045/i?3. Table SI summarizes the energies obtained by extrapolating 
the energies obtained using the second-order propagator for 3-4 different time steps to the zero time step limit. 
Comparison of Table II of the main text and Table SI shows that the energies for the model 2bZR-|-3bRp with 
p = 5 — 7 obtained by the two different approaches agree within error bars (for p = 4 and 8, we did not perform 
calculations using the fourth-order propagator). 

TABLE SI: PIMC energies for the model 2bZR+3bRp for = 5 — 15. Columns 2-6 show the scaled energy En/N/^Es/S) 
for p = 4-8, respectively. The error bars (not explicitly reported) are around 6%. 

N 2bZR-h3bR4 2bZR-h3bR5 2bZR-h3bR6 2bZR-h3bR7 2bZR+3bR8 


5 

5.58 

6.36 

6.67 

6.70 

6.63 

6 

8.22 

8.90 

9.53 

10.2 

10.4 

7 

8.96 

11.1 

12.1 

12.3 

13.4 

8 

10.4 

12.9 

14.8 

15.6 

16.2 

9 

11.6 

14.8 

16.8 

18.1 

18.7 

10 

12.7 

16.5 

17.9 

20.0 

21.2 

11 

13.1 

18.0 

20.5 

22.0 

23.2 

12 

13.5 

18.7 

22.0 

23.8 

25.3 

13 

13.5 

20.0 

23.5 

25.7 

27.6 

14 

14.3 

20.9 

24.8 

27.6 

29.1 

15 

14.6 

21.5 

25.9 

28.9 

30.5 


Table S2 reports our DMC energies, obtained by extrapolating energies for 4-5 finite time steps to the zero time 
step limit, for the models He-He(scale) and He-He(arctan) [2]. 

TABLE S2: DMC energies for the Hamiltonian with two-body van der Waals interactions for A = 4 — 15. Columns 2-3 show 
the scaled energy Em/N/[E s/Z) for the models He-He(scale) and He-He(arctan), respectively. The error bars (not explicitly 
reported) are around 1%. 

A He-He(scale) He-He(arctan) 


4 

3.952 

3.936 

5 

7.761 

7.737 

6 

11.81 

11.78 

7 

15.86 

15.80 

8 

19.76 

19.67 

9 

23.50 

23.38 

10 

27.06 

26.87 

11 

30.40 

30.11 

12 

33.60 

33.26 

13 

36.60 

36.18 

14 

39.45 

38.81 

15 

42.16 

41.66 
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